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Shape optimization of unconstrained and constrained damping 
layers is treated. The specific problem analyzed is a 
cantilever beam loaded at its tip by a harmonic force. 
Finite element modeling and mathematical programming 
techniques are used to obtain the solution. Performance 
measures are taken to be reduction of maximum displacement 
and increase in fatigue lifetime. Results include the 
improvement, over the uniform treatment case, of these 
measures when the profile of the damping layer is optimized. 


INTRODUCTION 

Treatment of vibration problems by damping layers, both constrained and uncon- 
strained is quite common. Early work in the field can be found in Ross, Ungar and 
Kerwin [1], More recently, finite element techniques have been used to address the 
problem. Papers relevant to the present work are those of Johnson, Kienholz and 
Rogers [2], Johnson and Kienholz [3], Soni and Bogner [4], and Soni [5]. 

Advances have also been made recently on the structural optimization front. 
Improvements in design sensitivity analysis were given by Kim, Anderson and 
Sandstrom [6]. Shape optimization techniques using pure finite element modeling 
were presented by Kikuchi, Chung, Torigaki and Taylor [7]. Of note is the work of 
Niordson [8], who showed the importance of imposing a slope constraint in the opti- 
mum design of elastic plates. Viscoelastic materials have also been treated. A 
study quite closely related to the theme of the present paper was given by Lekszycki 
and Olhoff [9], who analyzed shape optimization of an elastic beam covered by an 
unconstrained viscoelastic layer. Using calculus of variation techniques, they 
obtained an explicit optimality condition, which was solved in an iterative fashion. 

Here shape optimization is considered for both constrained and unconstrained 
layers on a beam. A key question addressed is, for a given volume of material, 
how much improvement can be obtained, over the uniform treatment case, if the profile 
of the damping layer is allowed to vary and be optimized. The specified problem 
treated involves a cantilever beam loaded at its tip by a time harmonic force. 
Performance measures are taken to be reduction in maximum displacement, and improve- 
ment in fatigue lifetime. Finite element modeling, together with numerical 
approaches to the complex eigenvalue problem and mathematical programming techniques 
are used to obtain the solution. 

It should be noted that complete details are not presented in the paper (these 
can be found in [10]). The work focuses on the essential ideas and on the results. 
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MECHANICAL MODELING 


Only a single constrained layer is treated here and the basic configuration 
is sketched in Fig. 1. The mechanical modeling used is traditional. The basic 



Figure 1. Basic Configuration of a Three-Layered Beam 


beam and constraining layer are taken to be Euler-Bernoulli beams and the damping 
layer, for which shear is important, is treated as a Timoshenko beam. Perfect 
bonding is assumed. With this modeling, the displacement components are given by 

U(x, z , t) = 


U 1 (x,t)-(z-d 1 )W x U,t), d 1 -T 1 /2<zld 1 +T 1 /2, base layer 
U 2 (x,t)-(z-d 2 )W x (x,t),d 2 -T 2 /2<z<d 2 +T 2 /2, covering layer 

U C (x,t)-(z-d C )l|> C (x,t), d C -H/2<z<d C +H/2, damping layer 
W (x, z, t) = W(x,t) 


(5) 


U C (x,t) = (l/2)(U 1 (x,t)-HJ 2 (x,t))+(l/4)(T 2 -T 1 )W 

ip C (x,t) = (l/H)(U 1 (x,t)-U 2 (x,t))-(l/2H)(T 1 +T 2 )W >x (6) 

In the above, U^, U 2> and U c are the midplane longitudinal displacements of the 
base beam, covering layer and core, respectively. 

The relevant, non-zero strains are 


xl 


= U. - (z-dJW 
l,x v 1 ,xx 


(7) 


£ x2 " U 2,x - < Z - d 2 )W ,xx 


( 8 ) 


< = (U l,x +U 2,x )/2+(T r T 2 )W ,xx /4 + (z - d c )[U 2,x- U l,x )/H + (T l + V W ,xx /2H] (9) 

( 10 ) 


y c = -^ C + w ,d C -H/2<z<d C +H/2 
xz * x 


For the base beam and covering layer, the stress strain relations are 
°xl = E i e xl 


a x2 E 2 £ x2 


(ID 

( 12 ) 


where E denotes Young’s modulus. The damping layer is treated as a Kelvin solid, 
for which the stress-strain relations are 


I 


X 




• c 

X 


(13) 


a 


c 


xz 


~c c 
G y 
'xz 



(14) 


where G° stands for the shear modulus of the core and r|^ C are material 

parameters characterizing the viscoelasticity. For harmonic loading, such as is 
being considered here, the complex modulus approach is adopted. Then 


c * C 

a c = e e 

(15) 

X X 

c * c 

0 = G Y 

xz xz 

(16) 


where 


* 

G = complex shear modulus of the damping layer 


= G c (i+in c ) 

(17) 
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E = complex Young's modulus of the damping layer 
= E C (l+i£ C ) 


(18) 


where n c and £ c are the loss factors for the damping material. For harmonic motion 
of frequency w, the relationships between the Kelvin parameters and the loss 
factors are 


G c = E C /(2(l+v C )) 


(19) 


n = w/g 


r = q w/e 1 


( 20 ) 

( 21 ) 


In the sequel, following Nashif, Jones and Henderson [11], Poisson's ratio v c 
is taken to be a constant. 

Equations (1) through (21) essentially set forth the mechanical modeling. 
The procedure then is straightforward. The principle of virtual work states: 


/(a <5e + a 6 v )dv + 6 V_ - 6 V = 0 

XX xz *xz I s 

V 


( 22 ) 


where 6Vj and 6v s denote virtual work by the inertia forces and surface tractions, 
respectively. Using eqs. (1) through (21) in eq. (22) leads to an integral expres- 
sion involving the "degrees of freedom" U^, U 2 , W, 3w/3x. This expression is then 
discretized using a finite element method. Note that similar modeling can be done 
for an arbitrary number of layers. 

FINITE ELEMENT MODELING 

Using eqs. (1) through ( 6 ), it can be shown that the volume integrals in 
eq. (22) reduce to line integrals in the x-direction. These integrals are then 
discretized using finite elements of length L^. Rod elements are used for axial 
displacements. Specifically, the shape functions are given by: 


U i (x) = [(l-x/L e ) x/L g ] [ (U 1 U ± 2 )] T 


(23) 


1 2 

where i = 1, 2 indicate base beam and covering layer, respectively, and IK , IK 
are nodal displacements. Beam elements are used to handle the transverse deforma- 
tions, with shape functions given by 

T 

e 2 ] (24) 


W (x) = [N x N 2 N 3 N 4 ][W j 


e 1 w 2 


„ _ 3W 

where 0 = 3 —, and 
ox 
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(25) 


2 2 3 2 

N = 1 - 3x /L + 2x / L 

1 e e 

N_ = x - 2 x 2 /L + x 3 /L 2 

2 e e 

N = 3x 2 /L 2 - 2x 3 /L 3 

3 e e 

2 3 2 

N. = - (x /L - x / L ) 

4 e e 

Again in eq. (24), superscripts indicate nodal quantities. 

Standard finite element methodology now applies. For harmonic forcing the 
procedure ultimately leads to, on assembling the various element matrices, 

- CO 2 [M] X + i[C] X + [K] X = F (26) 

where X is a vector of nodal parameters and F is a vector of nodal forces (magni- 
tudes). The stiffness, mass and damping matrices, [K], [M] , and [C] are lengthy, 
but straightforward expressions and will not be reproduced here. Note that the 
form i[C] for the damping matrix, which is frequency dependent, arises from use 
of the complex modulus approach. 


FATIGUE LIFE TIME CALCULATIONS 


Here the approach set forth in the SAE document. Ref. [12] is followed. In 
reality localized plastic flow occurs in fatigue and the nominal stresses and 
strains, a and e, should be replaced by the actual quantities S and e. Neuber 
introduced the following empirical rule 


e 


2 

a 


max 

SE 


(27) 


This equation has two unknowns, e and S, and the other needed relationship is 
the cyclic stress-strain curve for the material, which is curve-fitted by 


e 



(fr)"’ 


(28) 


where K' and n' are material parameters. Eqs. (27) and (28) are then solved 
iteratively to obtain S ma . The number of cycles to failure N f is calculated 
from another empirical ri¥ationship, namely: 


S 

max 


= a^(2N f ) b 


(29) 


t 

where and b are material parameters. Later in the paper an aluminum alloy 
(AL3015) is studied and for this material the parameters are 
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4 

E = 1.0 x 10 ksi 
K' =28.6 ksi 
n' = 0.093 
cr^ = 38.4 ksi 
b = -0.088 


OPTIMIZATION PROBLEM 

The optimization task is to find a vector b of design variables H^, i = 1,2.. 
...n, where H. is the thickness of the damping layer in the i^ finite element, 
which will minimize the objective function f, here taken to be: 

f : min (max | R ^ | ) j = 1 > 2 > N, 

where R. represents the deflection response at node j, subject to the constraints: 

volume constraint of damping layer - V = 0 

and inequality constraints: 

H. > 0 
x 

H. u - H i > 0, where h/ 1 is an upper bound for 
slope constraints: 

3h . 

H -I 77 —^ -I > 0, H a specified constant. (30) 

v 1 3x 1 v 

The constraint on the gradient in eq. (30) needs some explanation. The idea 
was introduced by Niordson [8] in a study on optimization of elastic plates. He 
pointed out that without it, exotic shapes (tending towards ribbed structures, 
with extremely thin stiffeners) would be generated. Apart from the practicality 
of such structures, the underlying theory (Kirchoff plate theory) is not valid 
for such rapidly varying shapes. To preserve the underlying theory, he restricted 
the design space to plates of slowing varying thickness, by means of a slope 
constraint. 

In the current work, the numerical approach requires that the constraints be 
differentiable. Hence the slope constraints are replaced by the equivalent state- 
ments : 
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(31) 


H 

v 


3H. 

+ — — > 0 
dx 


H - 

v 


3h. 

x 

IT 


> o 


(32) 


One remark should be made. For some unconstrained layers, it was found that 
the slope constraints were not necessary to obtain smooth shapes. However, the 
constraints were found to be essential for constrained layers. 

Mathematical programming techniques are used here to obtain the solution. A 
proven, reliable technique is employed. Belegundu and Arora [13] showed that the 
program SUMT has these features. Moreover, a listing is available in Kuester [14]. 

To use the program, sensitivity derivatives U^(b Q ) are required, where 


W 



(33) 


Partial differentiation of eq. (26) can be shown to give 


{- 0 ) 2 [M] + i [c] + [K]>U. = R 

“1 -C 


where 


R = ,n 2 3[M] x i ^ y 9 f K l v 

~ w 3bT - " 1 9bT - " W~ £ 


(34) 


(35) 


Note that in the present problems, [M] , [K] , and [C] have known analytical forms 
and the derivatives in eq. (35) can be carried out explicitly. Then eq. (34) has 
the same structure as eq. (26) and can be solved in the same fashion once the 
latter has been solved. 

The constraint equations in problem are simple algebraic expressions and their 
sensitivity derivatives can also be readily obtained. 


NUMERICAL STRATEGY 

Eigenvalue extraction was done by means of a subspace iteration technique 
together with Jacobi's method for matrix diagonalization. Response was calculated 
using a Gaussian direct elimination method, modified for complex equations. Design 
sensitivity coefficients were calculated as discussed in connection with eq. (35). 
Their magnitudes are then fed into the optimization scheme SUMT. 


PROGRAM VALIDATION AND RESULTS 

To check the accuracy of the finite element modeling, several calculations 
were done to determine the natural frequencies and loss factors and compared with 
results of Soni [5]. The comparisons involved a cantilever aluminum beam (7 inches 
long, .5 inches wide and .06 inches deep) with .06 inches thick aluminum face sheets. 
The material constants used were E = 1.0 x 10 ^ psi and p = 0.1 lb/in^. The core 
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c c 

material was ISD468. For this material, the behavior of G and X] with 
frequency is known. Poisson’s ratio V c was taken to have the constant value 
0.35. It was assumed (also assumed throughout the paper) that the loss factor 
is the same in dilation and shear, so that £ c = 0°. Tables 1 and 2 giv^^om- 
parisons of the first six undamped natural frequencies and the ratio f) 2 /n * 
r) 2 ^ r ' being the modal loss factor, respectively. Overall, quite good agreement 
is seen, lending confidence to the numerical procedures. 


Table 1. Comparison of Natural Frequencies 


Mode 

Number 

M. L. Soni [5] 

(hz) 

Present Result (hz) 


i 

64.70 


64.13 


2 

298.00 


296.80 


3 

748.20 


745.80 


4 

1409.50 


1403.70 


5 

2305.00 


2296.00 


6 

3447.00 


3400.00 


Table 

2. Comparison of 

the 

( r ) c 

Ratio n 2 /n 

Mode 

Number 

M. L. Soni [5] 


Present Result 


i 

0.2725 


0.2840 


2 

0.2401 


0.2450 


3 

0.1531 


0.1560 


4 

0.0878 


0.0896 


5 

0.0560 


0.0572 


The optimization phase of the program was checked on the following test 
problem given by Rosenbrock (see [14]): minimize the objective function f, where 

f = subject to 

constraints: 

0<x.<42 

0 < (x^+2x 2 +2x^)172 


Rosenbrock gave the solution: f = *3456.0, x-^ - 24, x 2 - 12, - 12. The present 

method gave f = -3453.8, = 23.4, x 2 * 12.2, x 3 * 12. Very good agreement is 

seen. This, and the fact that the trends in Lekszycki and Olhoff’s [9] work were 
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reproduced accurately (as will be seen shortly) led to the conclusion that the 
program was accurate. 

Results for unconstrained layers will now be given. 

The first material studied is the one used in [91, for which the parameters 
are: E c = 0.1x10^ lb/ in , r) c = 0.5, p c = 0.035 lb/ in and for the base layer 
= E c , p = 0.1 lb/in3. The dimensions are: T 1 = 0.06 inches, B (width) = 0.5 
inches, L = 7 inches (Soni f s example). 

Note that this damping material has quite large moduli (perhaps unrealisti- 
cally so). Moreover, note that, as in Ref. [9], the effects of shear are 
neglected. 

A harmonic force is applied at the tip. The thickness constraint H. u is 
taken to be 0.32 inches. The initial amount of damping material must be 
specified. A percentage measure is used, namely: 

a, i ^ volume of damping material 

% volume of damping material = = — 

volume of base layer 

Figure 2 shows a result for 100% damping using the present methodology, but 
ignoring shear effects. A symmetric configuration is used with equal amounts of 
damping material on the top and bottom of the beam. Only the upper layer is shown 
in the figure. In fact without a constraining layer, shear effects in the core 
are quite small and can always be neglected so that in effect an Euler- Bernoulli 
beam is used. The first bending frequency of the composite beam is (^ = 595 rad/ 
sec and the excitation frequency is a) = 20 rad/sec, so that we have a case of low 



Figure 2. Optimal Shape for High Modulus Material 
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frequency excitation. The same trend as in Ref. [9] is seen. It is interesting 
that the present results were obtained without having to specify a constraint on 
the slope (true for all the results on the unconstrained layers) . It was also 
found that the thickness constraint (0. 32 inches) was not active. 

Table 3 shows the improvements that can be obtained for various damping 
treatments. Two kinds of percent reduction in responses are defined by: 


RRU 


Max, response of bare beam - max, response with uniform damping 
Max. response of bare beam 


RRO 


Max, response of bare beam — max, response with optimal damping 
Max. response of bare beam 


Table 3. Performance Improvements 


% Volume 

RRU % 

RRO % 

Fatigue Life 
Uniform 

Fatigue Life 
Optimal 

100.0 

88.6 

93.0 

0.1775 x 10 5 

0.8764 x 10 5 

66.6 

79.9 

87.4 

0.6722 x 10 4 

0.4270 x 10 5 

33.3 

59.4 

72.7 

0.235 x 10 4 

0.909 x 10 4 

16.6 

38.1 

58.3 

0.129 x 10 4 

0.5994 x 10 4 


Optimization would seem to be worth the trouble. For example, for 33.3/4 damping 
an improvement (RRO-RRU) of 13.3% is seen. For 100% damping the improvement is 
4.4%. Better fatigue performance is also seen. Optimization led to 39.4%(log 
value) improvement for 33% damping and 74% improvement for 100% damping. 

High frequency excitation was also studied (not treated in Ref. [9]). 

Figure 3 shows the optimum shape for 0) = 750 rad/sec and 100% damping. It is 
interesting to note that the optimum profile has the opposite trend to that for 
the low frequency profile. 

As a next step in the study, a more realistic material was chosen, namely 
LORD-400. This is a medium shear modulus material with parameters obtained from 
Ref. [11]. The natural frequency now is = 220 rad/sec (based on 100% damping). 

The optimal shape for a low frequency excitation (a) = 20 rad/sec) is shown 
in Figure 4, for a 100% damping material and a bound on thickness of 0.24 inches. 
Note that the same shape trend is seen, as for the high shear modulus. The thick- 
ness changes in Figure 4 are severe and one may question the use of Euler-Bernoulli 
beam theory. A smaller upper bound on the thickness constraint was used, namely 
0.15 inches. The result is shown in Figure 5. A different, smoother shape is seen. 
This shape dependence on the thickness constraint was not observed for the high 
modulus material. 


404 




Thickness distribution 
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Figure 5. Optimal Shape (LORD-400) for Small Upper Bound 


Using the lower upper bound, the following values were found for performance 
improvement: at 100% damping, RRU = 9.0%, RRO = 19.0% ; at 66.6% damping RRU — 
4.7%, RRO = 12.1%; at 33.3% damping, RRU = 2.4%, RRO = 7.0%. Results on fatigue 
are: at 33.3% and 100% damping, the improvements are 5% and 12%, respectively. 

Though the gains are not as large as for the high modulus material, optimization 
still seems attractive. 

A typical result for high frequency excitation (u) = 240 rad/sec) is shown 
in Figure 6, for a thickness upper bound of 0.15 inches and 33% damping. 

Note that the same shape reversal as was seen for the high modulus material 
is found. RRO has the value 24.5% so optimization is also worthwhile at high 
frequencies. 

Constrained layers will now be treated (for LORD— 400) • The first item that 
should be mentioned is that now a slope constraint is required. . Figure 7 shows 
a shape obtained without such a constraint. Large oscillations in the profile 
are seen (as was seen by Niordson [8] in his work on elastic plates). Such shapes 
are not acceptable within the framework of the current mechanical modeling. It 
was discovered, like Niordson that a slope constraint led to smoother profiles. 

A slope constraint was imposed at the element level in the form 


H i + i - H i 


< .25 H 


i = 1,2, N-l 


where H is the thickness of the original uniform damping layer. The thickness 
of the covering layer is taken to be 10% of the thickness of the damping layer. 
Only symmetric configurations are considered. A thickness bound of .15 inches 
was used. 
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T hickness 



Figure 7. Optimum Shape With No Slope Constraint 
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An optimum shape for low frequency excitation (to = 20 rad/sec), is shown 

in Figure 8 for 100% damping. Comparing this with Figure 5, it is seen that 

the optimum profiles have opposite trends. One should not anticipate the same 

trend in both cases. The basic stress at work in the unconstrained damping layer 

is the bending stress a c , whereas it is the shear stress a c in the constrained 

x xz 

case. 

The improvement in performance (RRO) at 100% damping was found to be 53%. 

This should be compared with the 19% improvement noted for the unconstrained layer 
(with the smaller thickness bound). It can be concluded that constrained layers 
lead to significant improvement in performance. 


Line Chart for columns: X-j Y i ... X1Y2 
optimal □ uniform 



Figure 8. Optimum Shape of a Constrained Damping Layer 
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